% indexJ(3,4,[4,2,1])
% 
%      0     0     0
%      0     0     1
%      0     1     0
%      0     1     1
%      0     2     0
%      0     2     1
%      1     0     0
%      1     0     1
%      1     1     0
%      1     1     1
%      1     2     0
%      1     2     1
%      2     0     0
%      2     0     1
%      2     1     0
%      2     1     1
%      2     2     0
%      3     0     0
%      3     0     1
%      3     1     0
%      4     0     0
function ydot=testrdforce2(t,y,R,sig)

load('orthbase');

ydot(1) =R(t)-1/2*sig(t).^2;
ydot(2) =sig(t).*ms{3}(t);
ydot(3) =sig(t).*ms{2}(t);
ydot(4) =sig(t).*(ms{2}(t)+ms{3}(t));
ydot(5) =0;
ydot(6) =sig(t).*ms{3}(t);
ydot(7) =sig(t).*ms{1}(t);
ydot(8) =sig(t).*(ms{1}(t)+ms{3}(t));
ydot(9) =sig(t).*(ms{1}(t)+ms{2}(t));
ydot(10)=sig(t).*(ms{1}(t)+ms{2}(t)+ms{3}(t));
ydot(11)=sig(t).*ms{1}(t);
ydot(12)=sig(t).*(ms{1}(t)+ms{3}(t));
ydot(13)=0;
ydot(14)=sig(t).*ms{3}(t);
ydot(15)=sig(t).*ms{2}(t);
ydot(16)=sig(t).*(ms{2}(t)+ms{3}(t));
ydot(17)=0;
ydot(18)=0;
ydot(19)=sig(t).*ms{3}(t);
ydot(20)=sig(t).*+ms{2}(t);
ydot(21)=0;
ydot=ydot';